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ABSTRACT 

We use a Cartesian grid to simulate the flow of gas in a barred Galactic potential and inves¬ 
tigate the effects of varying the sound speed in the gas and the resolution of the grid. For all 
sound speeds and resolutions, streamlines closely follow closed orbits at large and small radii. 
At intermediate radii shocks arise and the streamlines shift between two families of closed or¬ 
bits. The point at which the shocks appear and the streamlines shift between orbit families 
depends strongly on sound speed and resolution. For sufficiently large values of these two 
parameters, the transfer happens at the cusped orbit as hypothesised by Binney et al. over two 
decades ago. For sufficiently high resolutions the flow downstream of the shocks becomes 
unsteady. If this unsteadiness is physical, as appears to be the case, it provides a promising 
explanation for the asymmetry in the observed distribution of CO. 
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1 INTRODUCTION 

More than twenty years ago, |Binney et al.| < [1991]) (hereafter 
BGSBU) proposed a consistent picture to explain spectral line 
emission by different species, HI, CO and CS, in the Galactic- 
centre region |/| < 10° and \b\ <2°. They related the flow of gas in 
an externally imposed, rigidly rotating barred potential to the struc¬ 
ture of the longitude-velocity, (/,v), plane that one obtains by pro¬ 
jecting closed orbits along lines of sight through the disc. Accord¬ 
ing to BGSBU, gas in the outer parts of the bar drifts slowly towards 
the centre while following x\ orbits. These orbits become more 
and more elongated as the centre is approached, and eventually the 
“cusped orbit” is reached, which has a cusp at each end. Interior 
to the cusped orbit, the orbits of the x\ family are self-intersecting. 
BGSBU hypothesised that when gas reaches the cusped orbit, it en¬ 
counters a shock, and then quickly plunges onto X 2 orbits. There¬ 
after, the gas drifts towards the centre following X 2 orbits. Using 
this simple representation of the gas flow, BGSBU provided an 
appealing interpretation of observational data: HI emission comes 
mainly from gas on non self-intersecting x\ orbits, CO forms as 
gas is shocked on reaching the cusped orbit, which explained the 
characteristic parallelogram-shaped envelope of CO emission in the 
(/,v) plane, and CS emission comes from dense, post-shocked gas 
flowing on X 2 orbits. 

The BGSBU picture has enjoyed considerable success and, 
alongside the photometric study of |Blitz & Spergel] j!991| >, con¬ 
vinced the community that our Galaxy is barred. In the following 
years, the presence of the bar has been confirmed by further photo¬ 
metric evidence ( [Stanek et al.|1994|[Dwek et al. |1 995 1 [Binney et al.| 
|1997| ), and there is now little doubt that the Milky Way is indeed 
barred. 

The BGSBU picture relied on the assumption that gas stream¬ 
lines follow closed orbits and that the x\ -A X 2 transfer happens 


at the cusped orbit, and needed validation by hydro simulations. In 
support of their model, BGSBU pointed to simulations by |Athanas^| 
|soula||T992b| ). However, these simulations used a different poten¬ 
tial from BGSBU, so a natural next step was to run hydro simula¬ 
tions in the potential BGSBU had used. One such study appeared 
( [Jenkins & Binney|1994) , but its results were not encouraging: in 
this simulation, orbits near the cusped x\ orbit, which played a cru¬ 
cial role in BGSBU picture, were found to be unoccupied by gas. 


More recently, numerous hydrodynamical simulations have 
been run with the goal of understanding the kinematics and dynam¬ 
ics of the cold gas in our Galaxy ([Mu lder & Liem||1986[ |Weiner 


& Sellwood||1999| |Englmaier & Gerhard||1999 |Lee et al"[fl 999 


Fux||1999[~ Bissantz et al.||2003[ |Rodriguez-Fernandez & Com bes 


2008[|Baba et al.|2010[|Pettitt et al.|2014| ). However, none of these 

simulations use the potential BGSBU used. Some used a poten¬ 
tial inferred from infrared photometric data (Englmaier & Gerhard 


|1999[|BTssantz et al.|2003[|Rodriguez-Fernandez & Combes|2 008), 

while others ( |Fux|1999| > used the potential generated by a combined 
hydro and N-body simulation of the Milky Way. Notwithstanding 
these efforts, the interpretation of the observational data remains 
problematic in some aspects ( [Sormani & Magorrian|2015j ). 


Most authors of papers using hydro simulations mention 
closed orbits, in a more or less explicit connection with BGSBU, 
but the literature lacks a detailed examination of the extent to which 
good hydro simulations support the BGSBU picture. The nearest 
the literature comes to filling this need is the paper of [Jenkins &[ 
Binney (1994), which describes simulations that are of low res¬ 
olution by today’s standard, and uses sticky particles rather than a 
conventional hydro simulation based on the Euler equations. Papers 
comparing closed orbits with the results of hydro simulations in bar 
potentials can be found |van Albada & Sanders) 1982[|Athanassoula| 
|1992a[b| , but, for lack of resolution or other reasons, none of them 
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provides sufficient detail to show which orbits are or are not occu¬ 
pied by gas, especially in the vicinity of the cusped orbit, which 
plays a crucial role in the BGSBU picture. 

In this work, we use high-resolution hydrodynamical simula¬ 
tions to re-examine the BGSBU picture. In the first part, we test 
the extent to which the physics of the gas flow hypothesised by 
BGSBU is supported by the simulations. Is the gas flow far from 
the shocks well approximated by closed orbits? Can the gas flow 
be understood as a transfer from x\ to X 2 orbits? Does the tran¬ 
sition happen at the cusped orbit as conjectured by BGSBU? We 
show that the answers these questions depend on the spatial reso¬ 
lution and sound speed used in a hydro simulation. The results are 
likely to be valid for all barred potentials that have a general resem¬ 
blance to the BGSBU potential. In the second part of the paper, we 
discuss the implications of our results for the interpretation of the 
observational data. Can we identify in the simulations structures 
reminiscent of the CO parallelogram of BGSBU? Under what con¬ 
ditions does the size of the X 2 disc match the region covered by CS 
emission? Can we explain the high velocity peaks in the HI (/,v) 
diagrams at |v| ~ 270kms -1 and |/| ~ 2°? 

This paper is structured as follows. In Sect. [2] we present the 
numerical schemes employed in the simulations. In Sect. [ 5 ] we show 
the results of the hydro simulations. We discuss the physical inter¬ 
pretation of the simulations in Sect. [4] and in Sect. [5] we discuss 
their implications for the interpretation of observational data. We 
finally summarise our findings in Sect. [6] 


2 METHODS 

2.1 Hydro Simulation Scheme 


We assume that the gas is a fluid governed by the Euler equations 
complemented by the equation of state of an isothermal ideal gas. 
Then we run two-dimensional hydrodynamical simulations in an 
externally imposed, rigidly rotating barred potential. The output of 
each simulation consists in snapshots of the velocity and surface 
density distributions p(x) and v(x) at chosen times. 

We use a grid-based, Eulerian code based on the second-order 
flux-splitting scheme developed by |van Albada et ah] (1982) and 
later used by |Athanassoula| ( [T992b| ), |Weiner & Sellwood| ( |1999| and 
others to study gas dynamics in bar potentials. We used the same 
implementation of the code as was used by |Sormani & Magor^j 
rian|(2015]>, s lightly modified to implement the recycling law of 


Athanassoula| ( |1992b| . 


This recycling law introduces a term in the continuity equation 
to take into account in a simple way the effects of star formation and 
stellar mass loss. The equation governing this process is: 


^ = a(pg - p 2 ) 


( 1 ) 


where a = 0.3 Mq pc -2 Gyr -1 is a constant and po is the initial sur¬ 
face density, which is taken to be po = lM 0 pc -2 . In practice, the 
only effect of the recycling law is to prevent too much gas accu¬ 
mulating in the very centre, and it does not affect the morphology 
of the results. Hence, the results of this paper do not change if we 
disable the recycling law. 

We used a grid N x N to simulate a square 10 kpc on a side. N 
depends on the resolution of the simulation. For example, if the grid 
cells are dx = 5 pc on a side, we have N = 2000. In each run the ini¬ 
tial conditions are as follows. We start with gas in equilibrium on 
circular orbits in an axisymmetrized bar and, to avoid transients, 
turn on the non-axisymmetric part of the potential gradually during 


the first 150Myr, in such a way that the total mass of the underly¬ 
ing potential is conserved in the process. We use outflow boundary 
conditions: gas can freely escape the simulated region, after which 
it is lost forever. The potential well is sufficiently deep, however, 
that very little gas escapes the regions of interest. 


2.2 The Potential 


We use the same potential as Jenkins & Binney (1994). This arises 
from two components. The first is the bar used by BGSBU, which 
has the density distribution: 


Pb(fl) = pbO 


(a/ao)~ a 

(a/a 0 ) _|3 


if a ^ ao 
if a > ao 


( 2 ) 


where a = yj'x 2 + (y 2 + z 2 )/q 2 and the values of the parameters are 
ao = 1.2kpc, a = 1.75, (3 = 3.5, Pbo = O.69M 0 pc -3 , q = 0.75, so 
the major axis of the bar always lies along the x axis. The sec¬ 
ond component is a razor-thin exponential disc, which has been 
added to complete realistically the circular velocity curve outside 
R ~ 1 kpc, and it has little influence inside this radius. The expo¬ 
nential disc is generated by a surface density distribution 

£(/?)= Loc~ r / r *, (3) 


where R is the radius in cylindrical coordinates and the parameters 
have values Eq = 1300M 0 pc -2 , Rj = 4.5kpc. 

The potential is assumed to be rigidly rotating with constant 
pattern speed D p = 63 km s -1 kpc -1 . This places the Inner Lind- 
blad Resonance at R \lr = 0.6kpc and corotation at Rq r = 3.7kpc. 


2.3 Projecting to the (/,v) plane 

We adopt a very simple projection procedure to produce the pre¬ 
dicted (/, v) distributions for each simulation snapshot (p(x), v(x)). 
Throughout this paper, we assume that the Sun is undergoing cir¬ 
cular motion at a radius Ro = 8 kpc with speed v 0 = 220km s -1 . 
Calling <|) the angle between the major axis and the Sun-GC line, 
the Cartesian coordinates of the Sun are given by x 0 = Rocosty, 
y 0 = Rosing- All the projections made in this paper assume <|) = 
20°. 

The resolution of our (/,v) diagrams is A/ = 0.25° in longi¬ 
tude and Av = 2.5 km s -1 in velocity. Along each line of sight, 
we sample the density and the velocity by linearly interpolating 
the results of the simulations at points separated by 5s = lpc. 
These density measures are accumulated in velocity bins of width 
Av = 2.5kms -1 . The final (/, v) intensity at the chosen longitude in 
each range of velocity is the mass in the relevant bin divided by the 
square of its distance. 

This procedure yields a predicted brightness temperature that 
is linear in column density so it is equivalent to the simplest ra¬ 
diative transfer calculation. In the case of HI, the brightness tem¬ 
perature is linear in the column density if the gas has constant 
spin temperature and its optical depth is negligible. So our pro¬ 
jection is equivalent to simple HI radiative transfer in the constant- 
temperature, optically-thin case. The assumption of constant tem¬ 
perature is known to be a simplification for Galactic HI, which is 
instead often modelled as a medium made by two or more phases at 
different temperatures (see for example [Ferriere|2001| ). In the case 
of 12 CO, the brightness temperature is not linearly related to den¬ 
sity when considering a single cloud, but a linear relationship will 
hold between brightness temperature and the number density of un¬ 
resolved CO clouds provided the cloud density is low enough for 
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Figure 1 . xi and X 2 orbits in the potential used by BGSBU. Orbits are shown 
in the frame that rotates (clockwise) with the bar. The horizontally elongated 
orbits form the x\ family, while the vertically elongated ones are from the 
X 2 family. BGSBU hypothesised that the orbits drawn in full lines have gas, 
while orbits those shown dashed are unoccupied. The cusped x\ orbit is the 
smallest horizontally elongated orbit drawn in full lines and is shown black. 


shadowing of clouds to be unimportant (see, e.g.,|Binney & Merri-| 
|field| 1998[ §8.1.4). 


3 RESULTS 

Fig.[l]shows a selection of closed orbits in the BGSBU potential, in 
the frame corotating with the bar. Orbits that BGSBU believed to 
carry gas are shown by full lines, while those they thought empty 
are shown by dashed lines. In the outer region we show a nested 
sequence of non self-intersecting x\ orbits that terminates in the 
cusped orbit, drawn in black. Inside this orbit we show dotted two 
self-intersecting x\ orbits. The vertically elongated orbits belong 
to the X 2 family, which extends quite a bit beyond the point where 
the cusped x\ orbit intercepts the vertical axis. BGSBU argued that 
gas transfers between the x\ and x 2 families at the cusped orbit. A 
primary goal of this paper is to test this conjecture. 

Fig. HI summarises the results of our simulations. It shows 
the density of hydro simulations for different grid spacings dx 
and sound speeds c s . All snapshots are taken at the same time 
t = 280Myr. Some common features of the gas flow can be iden¬ 
tified in all panels. In the outer part, approximately corresponding 
to the green region, the gas follows x\ orbits. At some point near 
the v axis, two thin offset shocks emerge, which connect the green 
region to the reddish central disc. The central disc is called the X 2 
disc and is made by gas on X 2 orbits. Most of the gas plunges to the 
X 2 disc or “central molecular zone” through the shocks. 

Thus, the gas follows the x\ orbits in the outer part and the 
X 2 orbits in the central part, with a transition zone containing the 
shocks in between. In each simulation, we can identify an inner¬ 
most occupied x\ orbit. The shocks are formed just after this orbit 
and they induce the transition from the x\ to the X 2 family. We call 
this innermost occupied x\ orbit the transition orbit, and the tran- 
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sition point its position in a parametrisation of the sequence of x\ 
orbits. 

BGSBU assumed that the transition orbit is the cusped orbit. 
In our simulations, the transition point and the size of the X 2 disc 
depend strongly on both the resolution and the sound speed. Con¬ 
sider, for example, the middle column in Fig. [2] corresponding to 
c s = lOkms -1 . As we increase the resolution, the transition point 
moves inwards while the X 2 disc shrinks. At the highest resolution, 
dx = 5 pc, the transition orbit almost coincides with the cusped or¬ 
bit, as predicted by BGSBU. At lower resolution the transition hap¬ 
pens earlier and the transition orbit is much bigger than the cusped 
orbit. Fig. [3] shows the same density snapshots as Fig. [2] superim¬ 
posed on the closed orbits that BGSBU thought carried gas for a 
better comparison. 

Increasing the sound speed also has the effect of postponing 
the transition and shrinking the X 2 disc. Consider, for example, the 
second row in Fig. [2] corresponding to dx = 20pc. At this resolu¬ 
tion, the transition happens very early, for c s = 5kms -1 , and the 
transition orbit is quite an outer x\ orbit (approximately the yel¬ 
low orbit in Fig. |TJ. As we increase the sound speed, the transition 
orbit moves inwards, and for c s = 20kms -1 it coincides in Fig. [I] 
with the green orbit that lies just outside the cusped orbit. We will 
discuss the origins of these systematics in Sect. [4] 

The two highest-resolution simulations in the right column 
of Fig. [ 2 ] look peculiar. These correspond to c s = 20kms -1 and 
dx = 5,10pc. At the time shown, all other simulations have already 
reached an approximate steady state, and they would not appear 
significantly different after another 250Myr or more. These two 
simulations instead manifest a complex unsteady flow interior to 
the transition orbit. The shocks are not stable, but keep forming 
and dissolving in an endless cycle. In some snapshots, they are al¬ 
most completely formed and smooth (see also Sect. |4.3| and Fig. [9j. 
At these moments, they lie very close to the x axis. A little later, 
vortices grow on the leading side, which move around and eventu¬ 
ally bump on the opposite shock, creating more vorticity. Later the 
cycle repeats. Tests described in Sect. [43] suggest that the unsteadi¬ 
ness is real and not an artifact of our particular code. 

So far we have not discussed the velocity structure of our sim¬ 
ulations, but this is crucial for the interpretation of observations. 
Fig. [i] shows the projections into the (/,v) plane of the snapshots 
of Fig. [2] for an assumed angle <|) = 20° between the bar’s major 
axis and the Sun-Galactic centre line. The lines in Fig. |4] show the 
(/,v) traces of some of the closed orbits plotted in Fig. fusing the 
same colour scheme - we plot only the orbits that BGSBU thought 
occupied (those shown with full lines in Fig. [TJ. 

For all resolutions and sound speeds, in Fig.[4]the envelope of 
the outermost x\ orbits matches the envelope of the hydro distribu¬ 
tion very well. As we move towards smaller values of |/|, the traces 
of orbits sweep up towards the high-velocity peak of the cusped or¬ 
bit, but at the transition orbit the hydro envelope starts to fall, and 
thus becomes separated from the orbit envelope. The projection of 
the hydro X 2 disc is clearly identifiable as a darker region near the 
centre. At the smallest longitudes its boundary is delineated by the 
traces of the X 2 orbits, but for low sound speeds or resolutions it 
extends far beyond the region covered by the plotted X 2 orbits. 

This finding again confirms the picture of the gas flow that 
we delineated above. When the transition point is too early, the 
innermost non-self intersecting x\ orbits, which are populated by 
gas in the BGSBU picture, are void of gas in the hydro simulations. 
Outer X 2 orbits that lack gas in the BGSBU picture are occupied by 
gas in the hydro simulation. As we increase the resolution at c s = 
lOkms -1 (middle column), the transition between the x\ and X 2 
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Figure 2. The fluid density in hydro simulations in the BGSBU potential for different spatial resolutions and sound speeds. c s is increasing left to right taking 
values 5,10,20kms _1 . dx is decreasing from top to bottom taking values 40,20,10,5pc. Gas has reached an approximately steady state in the rotating frame 
and circulates clockwise. All snapshots are taken at t = 280Myr. 

families moves inwards, and the hydro envelope matches more and envelope of the hydro the projection of the X 2 disc match very well 

more closely the predictions of the BGSBU picture. At the same as in the BGSBU picture, 

time, the projection of the X 2 disc shrinks, also approaching the 

BGSBU picture. Eventually, for dx = 5pc and c s = lOkms” 1 , the At low sound s P eed ( left column of Fig. [4), the transition al- 

ways happens early, so orbits close to the cusped orbit are unoc¬ 
cupied at all resolutions. Therefore at low sound speeds the hydro 
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Figure 3. Same as Fig.[2]but with superimposed orbits. 
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Figure 4. The simulations of Figs.[2]and[3]projected into the (/, v) plane. Solid lines show the (/, v) traces of closed orbits, with colours matching other figures. 
The Sun is assumed to be in a circular orbit with v = 220kms -1 at Ro = 8kpc, and the bar major axis makes an angle <|> = 20° with the Sun-Galactic centre 
line. 


simulations gives results inconsistent with the BGSBU picture. At 
high sound speeds the hydro simulations approach the BGSBU pic¬ 
ture as we increase the resolution at first. But at higher resolution, 
unsteadiness makes the projected hydro deviate significantly from 
the BGSBU picture. 


To explain this situation from a face-on perspective, we take 
as an illustrative example the case dx = 10, c s = lOkms -1 . Fig. [ 5 ] 
compares the velocity field of the hydro simulation with that of the 
closed orbits. The left panel shows the velocity field of the hydro 
simulation. The central panel shows the best approximation to the 
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x[kpc] x[kpc] x[kpc] 

Figure 5. Using closed orbits to approximate the hydro velocity field for the simulation with dx = lOpc, c s = lOkms -1 shown at left. The central panel shows 
the best approximation to that velocity field obtained using closed orbits belonging to x\ or X 2 family. For each point where more than one orbit is found, 
the chosen velocity is the one closest to the hydro velocity. The right panel shows the vector differences between the left and middle panels. In all panels 
points without orbits are masked and the density of the hydro simulation is visible in the background. The colorbars show the speed at each point in units of 
100 kms -1 . 


hydro velocity field that can be obtained from the x\ and X 2 or¬ 
bits: at each point where more than an orbit is present, we show 
the velocity of the orbit that best matches the hydro velocity field. 
In all panels the density of the hydro simulation is shown in the 
background. Locations through which no closed orbit passes are 
masked out. The right panel is the most interesting panel: it shows 
the vector difference between the left and middle panels, and shows 
clearly that the outer x\ orbits and the inner X 2 orbits both repro¬ 
duce the hydro velocity field accurately. The red orbit in Fig. [T] is 
the transition orbit in this case, the last orbit at which the velocity 
fields coincide, and is just outside the shocks. After this, the shocks 
emerge and the two velocity fields suddenly diverge. 

At different resolutions and sound speeds, the situation is qual¬ 
itatively very similar but the point of transition between orbit fami¬ 
lies changes. For c s = 20km s -1 and high resolution the hydro flow 
in the transition region is unsteady. 

4 THE PHYSICS OF THE GAS FLOW 
4.1 dependence on the sound speed 

In Sect. [3] we have seen that as the sound speed increases, the tran¬ 
sition point and the shocks move inwards. At the same time, the 
shocks also become more horizontal and closer to the v axis. En- 
Iglmaier & Gerhard| < |1997| ) and |Patsis & Athanassoula| ( [2000] ) found 
similar results when varying the sound speed. A qualitative expla¬ 
nation for this behaviour is as follows (for another discussion see 
also [Englmaier & Gerhard] 1991 ) . 

The more elongated an x\ orbit is, the larger the variation 
along the orbit in the speed of an orbiting particle. Since the speed 
decreases as particles move from the orbit’s minor axis to its major 
axis, the density of gas that is streaming along the orbit increases as 
the major axis is approached. When the sound speed is high, pres¬ 
sure assists gravity in slowing gas as it approaches the major axis, 
and equally assists gravity in accelerating the flow after the major 
axis has been passed. In the absence of a shock, the gas is reversibly 
compressed and decompressed so it can continue to keep close to 
one orbit for several revolutions. 

If the sound speed is too low, the convergence of the gas as 
the major axis is approached leads to shock formation. Entropy is 
created in the shock, so the decompression after the major axis has 


been passed does not reverse what happened as the axis was ap¬ 
proached, and the flow deviates strongly from orbits. 

The key to avoiding shock formation is the ability of sound 
waves to carry information about fluctuations in density upstream 
so oncoming gas can be slowed in a timely manner when the den¬ 
sity increases at a downstream location. The nearer the cusped orbit 
is approached, the larger is the velocity gradient up which sound 
waves have to travel if a shock is to be avoided. Hence decreasing 
the sound speed causes the shock to form further out. 


4.2 Dependence on numerical resolution 

Increasing the resolution moves the orbit at which shocks form in¬ 
wards with important consequences for the interpretation of the gas 
flow in our Galaxy that will be discussed in Sect. [5] A likely expla¬ 
nation for this phenomenon is that at low resolution the innermost 
x\ orbits are inadequately resolved. To resolve the cusp we need in 
principle infinite resolution. By increasing the resolution, we can 
resolve orbits closer and closer to the cusped one, and gas can set¬ 
tle on these orbits, delaying the formation of shocks. 

Fig[6]tests this idea by showing the (/, v) traces of nine orbits 
from the cusped orbit (at the top) outwards. Each orbit is mapped 
four times at resolutions that increase from left to right. The faint 
solid lines show the traces of the orbits and as such are the same 
along every row. The big points scattered around this line show 
the (/,v) trace one obtains by associating each hydro cell through 
which the orbit passes with the mean velocity of its passage, and 
then projecting the resulting velocities of the visited cell. The small 
points show the projection at the visited cells of the hydro velocity 
field for c s = lOkms -1 . 

We see that orbits that are well outside the cusped orbit (lower 
part of the figure) are well represented by even the coarsest grid 
(left side of the figure), and, moreover, their velocities are accu¬ 
rately reproduced by the hydro code. As we move up the figure and 
therefore approach the cusped orbit, at low resolutions the big dots 
scatter widely around the orbit’s curve. The scatter is widest at the 
largest values of |/| because that is where the velocity gradient is 
largest and thus the finite resolution of the grid has its biggest im¬ 
pact. For each orbit outside the cusped orbit there is a resolution 
finer than which the hydro velocity field is almost the same as the 
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Figure 6. Each panel refers to a particular x\ orbit. In full lines, the trace of the orbit. The bigger, empty dots are the velocity field of closed x\ orbits sampled 
with the same resolution of a hydro simulation along the xy trajectory of the orbit. Smaller, filled dots are the projection of a hydro simulation cells along the 
same trajectory. Each column refers to one hydro simulations. From left to right, the resolution increases from dx = 40pc to dx = 5pc. All simulations are for 
c s = lOkms -1 . 
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orbital velocity because along that orbit the flow does not encounter 
a shock. 


4.3 Unsteady flow 

As mentioned in Sect. [3] the flow becomes unsteady at high sound 
speed and high resolution. To test whether this phenomenon was an 
artifact generated by our code, we re-ran some simulations with a 
completely different code, PLUTO ( [Mignone et al.|2007| . This is 
free software for the numerical solution of systems of conservation 
laws targeting high Mach number flows in astrophysical fluid dy¬ 
namics. It has a modular structure that makes it easy to change the 
algorithm used to simulate a flow. We use this modularity to inves¬ 
tigate the effects of the choice of Riemann solver and flux limiter. 
In all of our runs we used a static Cartesian mesh and RK2 time¬ 
stepping. The code was modified to implement our recycling law 

Q- 

We found that the Roe and HLL Riemann solvers produced 
the same results as each other, and our own code. The choice of 
flux limiter, by contrast, does have a significant impact on the com¬ 
puted flow: unsteadiness of the flow can be suppressed by changing 
to a more diffusive flux limiter. Flux limiters are used in numerical 
schemes to handle the flow close to discontinuities such as shocks. 
We tested three different limiters and found that unsteadiness is 
suppressed by choosing a more diffusive limiter, which has a higher 
numerical viscosity and produces thicker shocks. Fig. [7] shows the 
result of using PLUTO to perform a standard one-dimensional test 
problem, the SOD shock tube ( [Sod|1978| . It shows that the three 
limiters produce shocks of different thicknesses. The most diffu¬ 
sive (more viscosity) is the MINMOD limiter (Roe 1986 ). The Van 
Albada limiter used in this paper ( [van Albada et al.[1982| ) has inter¬ 
mediate diffusivity, while the least diffusive limiter (least viscosity) 
is the MC limiter ( [van Leer|1977| ). Fig. [8] shows a snapshot at an 
intermediate time for simulations with dx = 5 pc, c s = 20km s -1 
obtained using PLUTO. The only difference between the three sim¬ 
ulations is the limiter used. From top to bottom the figure shows the 
flows obtained with limiters of increasing diffusivity. The top flow 
has much irregular and unsteady structure that is completely ab¬ 
sent from the bottom flow. The central panel shows an intermediate 
level of unsteady structure. 

In these flows unsteady features seem to arise at the shocks, 
which develop wrinkles that shed vortices. The vortices then move 
away from the shocks. Vorticity is generated at the shocks and prop¬ 
agates away from them into the body of the flow. 

Fig.[9]shows the flow computed with the most diffusive limiter 
(MINMOD) at three different times. It shows the shocks cyclically 
straightening out and then developing wrinkles. More unstable sim¬ 
ulations also display a cyclical tendency for the shocks to straighten 
out and dissolve, but in these lower-diffusivity calculations the dis¬ 
solution of the shocks is more sudden. From these tests we con¬ 
clude that unsteadiness is not a feature of our particular hydro code, 
but is shared by other well tested codes as well. 

Turbulence is a consequence of high rates of shear in a high 
Reynolds number flow. In the flows studied by engineers, for ex¬ 
ample in pipes and over aircraft wings, turbulence arises in the thin 
boundary layer that forms when a low-viscosity fluid meets a solid 
surface. Vortices formed in the boundary layer move into the bulk 
of the fluid, making the flow generally unsteady. In our simulations, 
vortices form along the shocks (see |Binney|1974[ for an analysis of 
vorticity generation in shocks), where the shear rate is exception¬ 
ally high. The rate of shear is inversely proportional to the width 
of the shock, so it increases with the grid resolution and decreases 



Figure 7. Results of SOD shock tube test problem at t = 0.2 for different 
flux limiters. This figure zooms around a shock. The resolution is dx = 
0.0025. 


with the diffusivity of the limiter. Thus simulations with the finest 
grids and the least diffusive limiters have the highest shear rates 
and are most likely to become turbulent. 

The fact that turbulence appears only in simulations with a 
high sound speed can also be explained by the connection between 
shear rate and the onset of turbulence: the rate of shearing increases 
along the sequence of x\ orbits as one approaches the cusped orbit, 
so at low sound speeds, when the shock forms far from the cusped 
orbit, the maximum rate of shearing is smaller than when the sound 
speed is larger. 

A good indicator of the amount of shear in a 2D flow is the 
quantity (see for example |Maciejewski|2008| ) 


/ dv x dv y \ 

\dy dx ) \dx By) 


(4) 


This quantity is invariant under rotations of the coordinates, being 
the magnitude of the eigenvalues of the traceless shear tensor, de¬ 
fined by 


Dij ~ 2 (dx] ' dxi S,7(V ' V 


(5) 


Fig. [TO] shows the quantity xdx, where the shear x is estimated by 
finite differences. As claimed above, the shearing rate is high along 
the shocks and is higher for higher sound speeds. 

[Kim et alT|(2012) also encountered the onset of turbulence as 
isothermal gas moves in a rotating barred potential, and our Fig 
[8] is similar to their Fig.4, panel (c). They point out that vortices 
generated in one shock pass to the shock on the other side of the 
galaxy, and are there amplified. This process has been called the 
wiggle instability (see for example Wada & Koda 2004, Kim et al.| 
|2014[ and references therein). There is still no consensus on the 
nature of this instability: Wada & Koda (200 4]) conject ured it is a 
Kelvin-Helmholtz (KH) type instability, but Kim et al.^ < |2014| ) use 
a shearing box analysis to argue it should be considered as a differ¬ 
ent type of instability. The instability has also been attributed to a 
numerical noise caused by the discretisation of the fluid equations 
( jHanawa & Kikuchi|2012| ). 
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Figure 10. The quantity xdx, where x = is an indicator of the shear and is defined in Eq. {4}. In each panel, the maximum value reached by the quantity is 
shown. 
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Figure 8. Three snapshots capturing the moment when the instability is 
forming. The three snapshots come from three different simulations ob¬ 
tained with the Pluto code. The only difference between the three is the 
limiter used. On top using the MC limiter, middle the VanAlbada limiter 
and bottom the MidMod limiter. 
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Figure 9. Three late snapshots of a simulation using Pluto code and the 
MinMod limiter. This limiter allows the simulation to be only weakly un¬ 
stable. It shows that shocks are almost formed and destroyed cyclically. 
Swirls produced at one shock can propagate and bump on the shock on the 
other side. 
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Figure 11. The vertical radius of the X 2 disc as a function of resolution. Each 
point is obtained averaging the radius over 4 different snapshots at times 
t = 274,313,352,391 Myr. For resolution dx = 2.5pc the flow is unsteady 
and the size of the disc fluctuates, and the red crosses show the scatter in 
the size of the X 2 disc over the four snapshots. The horizontal dashed line 
shows the value of the intercept of the cusped orbit with the vertical axes. 


4.4 Numerical convergence 

We have seen that as we increase the resolution at given sound 
speed, the shocks are postponed and the %2 disc shrinks. Have we 
converged, or would a further increase in resolution produce sig¬ 
nificant changes in the flow? In order to discuss this, let us con¬ 
sider how the size of the x 2 disc depends on the resolution. In each 
snapshot in our simulations the vertical edges of the X 2 disc can be 
clearly identified by a sharp variation in the gas density as we move 
along the vertical axis. Consequently, the vertical size of the X 2 disc 
is well defined. 

Fig .[IT] shows the vertical radius of the X 2 disc plotted against 
resolution for simulations with sound speed c s = lOkms -1 . The 
data points for dx ^ 5 pc are from the simulations presented above, 
while the black dot at dx = 2.5 pc is obtained for an additional simu¬ 
lation. Since the computational cost of a simulation of given spatial 
extent rises as (dx) -3 , the additional simulation was run on a grid 
that covered an area only 5 kpc x 5 kpc in extent - we have verified 
that the results of lower resolution simulations for the region inside 
the shocks are unaffected by a reduction of the extent of the grid, so 
the additional simulation should provide a valid additional datum 
for the radius of the X 2 disc. 

The simulation with dx = 2.5 pc was unsteady as expected, 
given the tendency for unsteadiness to arise at high resolutions and 
sound speeds. On account of the unsteadiness, the size of the X 2 
disc in the dx = 2.5 pc simulation fluctuates. In Fig. |TT] the point 
for the dx = 2.5 pc simulation is obtained by averaging the X 2 disc 
size over four different snapshots at times t = 274,313,352,391 
Myr, and the red crosses show the values in each snapshot. Such 
averaging would not modify the data points for the coarser, steady, 
simulations. Since in the simulation with dx = 2.5 pc the mean ra¬ 
dius of the X 2 disc is compatible with this radius when dx = 5 pc, 
we conclude that the size of the X 2 disc has converged to near its 
true value at c s = lOkms -1 . 

As can be seen from Fig. |TT] the value of the vertical radius 
towards which we have converged for c s = lOkms -1 is the value 
of the intercept of the cusped orbit with the vertical axis. Inspection 
of Fig. [2] shows that for higher sound speeds, the X 2 disc becomes 
smaller. Indeed, [Englmaier & Gerhard] ( j !997| ) found that the X 2 disc 
disappears for high enough sound speed. 


4.5 A paradox 

In the limit of vanishing pressure, the characteristics of the Eu¬ 
ler equations reduce to ballistic trajectories. Hence, we expect 
the ballistic approximation to fluid flow to work best when the 
sound speed is low. Why is the BGSBU picture, which is founded 
on the ballistic approximation, more accurately reproduced with 
c s = lOkms -1 than with c s = 5kms -1 ? 

The answer is that the BGSBU picture includes the require¬ 
ment that the transition between orbit families lies close to the 
cusped orbit. This transition is effected by shocks, which are more 
rather than less likely to form in the limit of vanishing pressure. 
Shocks reflect the necessity in a fluid for there to be a unique ve¬ 
locity at each spatial point. As |Riemann] ( | 1 860| ) showed, in certain 
circumstances the Euler equations predict more than one velocity at 
a given point of a fluid in the sense that characteristics of the same 
family can cross. The unphysical implications of crossing charac¬ 
teristics are voided by a shock forming near the crossing point. In 
the shock the finite mean-free path of the fluid’s constituent parti¬ 
cles becomes important, and the Euler equations, which are based 
in the fiction that the fluid forms a continuum, cease to be valid. 

In a region of the disc where there are two distinct families of 
closed orbits, the principle that in the limit of vanishing pressure 
the streamlines/characteristics will be closed orbits does not suffice 
to determine the flow. One needs in addition some way of selecting 
one orbit family for the streamlines. Since the x\ orbits do not exist 
at very small radii and the X 2 orbits do not exist at very large radii, 
it is evident that in a barred galaxy the flow must at some point tran¬ 
sition from x\ to the X 2 orbits. Our simulations have shown that the 
location of the transition is controlled by the pressure even though 
the latter is too weak to be dynamically significant upstream of the 
transition. BGSBU argued that the transition occurs at the cusped 
orbit, the latest that it could occur, and with hindsight it is not sur¬ 
prising that this requires that the sound speed is > lOkms -1 . 


4.6 Sound speed and the nature of the ISM 

Our results can reproduce BGSBU predictions only for a particular 
sound speed. Does this mean we are providing a measurement of 
the sound speed? 

In addressing this question, a major issue is that we have in¬ 
vestigated only one potential and only a single value of the pattern 
speed, so we cannot exclude the possibility that in another potential 
the data could be fitted with a lower value of c s . A trivial illustra¬ 
tion of this fact is provided by one of the axisymmetric potentials 
that were used to interpret (/, v) plots before BGSBU. However, the 
sharply peaking circular-speed curves that these potentials require 
to fit the HI envelope in the (/, v) plane are highly implausible, and 
if the mechanism proposed by BGSBU for generating this feature 
is to work, gas must stay close to x\ orbits until rather close to 
the cusped orbit. That is, in a model with a small value of c s , the 
upward rise in the HI envelope would have to be at least partly 
explained by inward-increasing values of the circular speed. More¬ 
over, the cusped orbit would have to occur at smaller radii than 
BGSBU hypothesised. 

Given the almost fractal nature of the ISM, the physical sig¬ 
nificance of the parameter c s that plays such an important role in 
the simulations is far from evident. It is, however, worth noting that 
we have pushed the resolution so high that the cell size of our sim¬ 
ulations is comparable to the size of a giant molecular cloud - even 
small GMCs are > 5 pc across. 
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5 IMPLICATIONS FOR THE INTERPRETATION OF 
OBSERVATIONAL DATA 

5.1 BGSBU revisited 

The three key points in BGSBU’s interpretation of the observa¬ 
tional (Z,v) plots were: 

(i) The high velocities, reaching v ~ 270km s -1 , found in the 
HI (/,v) data could be explained by gas moving on x\ orbits just 
outside the cusped orbit. When the (/, v) plots had been interpreted 
in the context of an axisymmetric Galaxy in previous studies (see 
for example [Sofue|2013| ), the high velocities observed had required 
the circular-speed curve to first rise, and then fall with implausi¬ 
ble rapidity. By hypothesising a bar, BGSBU were able to explain 
the sharp peaks in the (/,v) envelope with a monotonically rising 
circular-speed curve. 

(ii) The parallelogram they identified in the CO emission was a 
cut version of the longitude-velocity trace of the cusped orbit. The 
discrepancy between the cusped orbit and the CO parallelogram 
was attributed to the limits of the ballistic approximation, and it 
was hoped that full hydrodynamical calculation would resolve this 
discrepancy. In particular, BGSBU pointed to the hydrodynamical 
simulations of Athanassoula (1992b ), which used the same numer¬ 
ical method as this paper, to construct an qualitative argument that 
qualitatively resolved the discrepancy. 

(iii) The emission at low longitudes visible in the CS was at¬ 
tributed to gas on X 2 orbits. Outer X 2 orbits that can be found in the 
potential (see Fig. [TJ were assumed to be unoccupied by gas. 

In the previous section we confirmed that the flow of gas can be un¬ 
derstood as a transfer from x\ to X 2 orbits. However, the transition 
point is not always at the cusped orbit, as BGSBU assumed, but de¬ 
pends on the parameters of the gas flow, especially the sound speed. 
For the right sound speed, c s = lOkms -1 , and our maximum res¬ 
olution, dx = 5 pc, the gas flow is very similar to that hypothesised 
by BGSBU. We now re-analyse the observational data in light of 
this particular simulation. Further below we take a step back and 
discuss interpretation of data in a broader view. 

In the centre-bottom panel of Fig. [4] we can see the (/, v) plot 
corresponding to the simulation with dx = 5pc and c s = lOkms -1 . 
The envelope of the hydro follows that of the closed orbits very 
well and reproduces the high velocity peaks found in observations. 
Fig. 1 12| shows how the line-of-sight velocity varies in the xy plane 
for the assumed observing angle <|> = 20°. The black dots show the 
locations that generate the envelope of emission in the (/,v) plane. 
We can see that the high-velocity peaks at |/| ~ 2° arise from gas 
that lies very nearly along an x\ orbit just outside the cusped orbit, 
as well as the portion of envelope that runs from the peak down 
into the zone of forbidden velocities. Thus item (i) above is nicely 
confirmed by this hydro simulation. 

Well outside the cusped orbit the black dots in Fig.[3]have wig¬ 
gles associated with spiral arms that run out from the extremities of 
the bar. These wiggles are reflected in small oscillations in the (/, v) 
envelope, qualitatively similar to what is observed for the envelope 
of real observations. 

As mentioned above, two distinct structures resembling a par¬ 
allelogram play a role in the BGSBU picture: one is the CO paral¬ 
lelogram, presumed to be a cut version of the cusped x\ orbit, and 
a second is the trace of an x\ orbit just outside the cusped orbit, 
which is responsible velocity peaks in the envelope of HI emis¬ 
sion. In the simulated (/,v) plots we can identify only one such 
parallelogram, namely the second one. How can we solve this in¬ 
consistency? Fig. [13] shows the shocks (blue and green) and the 



X[kpc] 

Figure 12. The distribution of the projected line-of-sight velocity in the xy 
plane for the model dx = 05pc, c s = lOkms -1 . The black dots show the 
points corresponding to the envelope of gas in the (Z,v) plane. The angle 
between the Sun-Galactic centre line and the bar major axis is assumed to 
be <|> = 20°. The colorbar is in units of lOOkms -1 . 

X 2 disc (red) in the xy plane (top) and in the (Z,v) plane (below). 
In the lower panel the shocks closely follow the vertical sides of 
the cusped orbit’s (black) parallelogram, but are cut almost exactly 
where the CO parallelogram ends. We conclude that the vertical 
sides of the CO parallelogram are formed by shocks. The other 
two sides of the CO parallelogram must be made up of gas flow¬ 
ing from one shock to the other. The shocks do not show up in our 
(/,v) projection of the gas flow, though they are quite apparent in 
the xy plane. In reality, we expect the shocked gas to be brighter 
than our simple minded radiative calculation suggests, because a 
lot of atomic gas is converts to molecular as it is compressed at 
the shocks. When this effect is taken into account the second paral¬ 
lelogram should appear in the (Z,v) projection of the model. Thus 
we confirm item (ii) above of BGSBU’s interpretation, although the 
physical mechanism that causes the parallelogram to be cut was not 
exactly described by BGSBU. Finally, Fig.[l3]clearly shows that in 
the (/, v) plane the central molecular zone (red) occupies the region 
of the inner X 2 orbits, as conjectured by BGSBU [item (iii) above]. 

When the x\/x 2 transition happens well before the cusped or¬ 
bit because either the sound speed or the resolution is low, there is 
no gas able to explain the high velocity peaks as in item (i). More¬ 
over, the shocks project to \l\ >2° so they cannot be associated 
with part of the CO parallelogram, and the central molecular zone 
is predicted to extend to higher \l\ than where significant CS emis¬ 
sion is found, compromising the interpretation of item (iii). Thus 
several different aspects of the data point to gas remaining on x\ 
orbits right up to the cusped orbit. 

5.2 The asymmetry 

It is well known that the molecular emission in the central molec¬ 
ular zone is highly asymmetric: three quarters of the 13 CO and CS 
emission comes from positive longitudes. Perspective effects can¬ 
not account for this asymmetry ( [Jenkins & Binney||1994| ), so the 
observed asymmetry must be transient: observations made tens of 
megayears in the past or future would often show asymmetry in 
the opposite sense. Thus the likely explanation of the asymmetry 
is unsteady flow, and the principal motivation of the work reported 
in [Jenkins & Binney|jl994j ) was to find evidence of unsteady flow. 
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Figure 13. Top: points on the narrow shocks and the X 2 disc in the xy plane. 
Bottom: their projected position in the (/, v) plane. In black the trace of the 
cusped orbit is shown. 

Although low-amplitude unsteadiness was generated in their sim¬ 
ulations, the present study implies that their simulations were too 
crude to probe the physically interesting regime. 

We saw above that the greater prominence of the CO rectangle 
in the observed (/, v) plot than in the theoretical (/, v) plot suggests 
that the shocks are important sites for the conversion of atomic to 
molecular gas. This being so, unsteady flow through the shocks will 
give rise to unsteady conversion of atomic to molecular gas, so the 
atomic/molecular ratio on each side of the Centre could well fluctu¬ 
ate as widely as the observations imply. However, a full explanation 
of the observed asymmetries in molecular-gas emission must await 
high-resolution simulations that keep track of the chemistry of the 
ISM. 

5.3 What do we still need to explain in the (/, v) diagram? 

Two aspects of observations remain inadequately explained: 

• Coherent broad features like the 3kpc arm and its counter¬ 
part on the far side of the Galaxy [Dame & Thaddeus[f2008] ). These 
are not produced in our simulations. They are produced in other 


hydro simulations |Mulder & Liem|l986l|Rodriguez-Femandez &| 
[Combes 2008), but these simulations did not reproduce the high 
peaks in the envelope of HI emission in the (Z,v) plane, probably 
for lack of resolution. 

• Forbidden emission at large longitudes. The portion of the 
(/,v) diagram covered by forbidden emission in our simulations 
is smaller than the region in which coherent forbidden emission is 
seen in the data. A higher quadrupole moment is probably needed 
to reproduce this. 

The essential elements of the BGSBU picture are that streamlines 
coincide with x\ and X 2 orbits, and that the shocks responsible for 
the transition lie near the cusped orbit. BGSBU illustrated these 
principles with one particular, very simple potential. Better fits to 
the data could surely be obtained with other, similar potentials. A 
fast way to select potentials worthy of closer examination would be 
to use closed orbits as BGSBU did. 


5.4 Relation to prior work 

The question of what physical mechanism determines the size of 
the X 2 disc is relevant for the interpretation of observations in our 
and external galaxies (see for example Combes 1996, Kim et alT] 
|2012| >. Our results suggest that some previous studies may be bi¬ 
ased by not taking into account the effects of varying the resolu¬ 
tion. For example, [Cole et al7| ( |2014| ) in their simulations of galaxy 
formation found that the main mismatch between their models and 
the observations was that the nuclear discs of their models were too 
big relative to their bars. Since nuclear discs are X 2 discs, our re¬ 
sults suggest that the mismatch would be resolved by an increase 
in resolution. 

Finally, we mention that our findings do support the hopes of 
|Bissantz et al.| ( |2003| . These authors found, as in our low-resolution 
simulations, that innermost non-self intersecting x\ orbits near the 
cusped orbit were unoccupied, and attributed this fact to details of 
the SPH scheme, which do not apply for our grid-based simula¬ 
tions. This work gives strong support to their view that in higher- 
resolution, grid-based simulations their inner x\ orbits would be 
occupied by gas, giving rise to peaks in their (/, v) projections. 


6 CONCLUSION 

|Binney et al.| ( [1991| (BGSBU) constructed a picture of the flow of 
gas through the central few kiloparsecs of our Galaxy. Their picture 
was based on the idea that gas follows closed orbits, and it involved 
a particular choice of orbit at which the gas transitions from the x\ 
orbit family to the X 2 family. Their orbit-based picture required val¬ 
idation by hydrodynamical simulations of gas flow. Early efforts in 
this direction did not provide the necessary validation, in part be¬ 
cause they did not adopt the same Galactic potential as BGSBU, but 
largely because in them gas occupied some X 2 orbits that BGSBU 
required to be empty, and left empty some x\ orbits that BGSBU 
required to be occupied. We have run high-resolution, grid based 
hydro simulations of gas flow in the potential of BGSBU and vali¬ 
dated their picture in the case that the effective sound speed in the 
ISM is c s > lOkms -1 . 

The simulations confirm that, regardless of the sound speed 
adopted and the grid resolution employed, gas streamlines closely 
coincide with closed orbits everywhere outside a shock-dominated 
transition region that divides the outer region, in which gas fol¬ 
lows x\ orbits, from an inner region, in which gas follows X 2 orbits. 
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However, the orbit at which the shock arises, and the transfer com¬ 
mences, depends on both the sound speed and the grid’s resolution. 
Shock formation is favoured by both low sound speeds and low 
grid resolutions, so increasing the sound speed and/or the grid res¬ 
olution moves inward the shock that causes gas to plunge from x\ to 
X 2 orbits. The BGSBU picture calls for the shock to occur as close 
to the Galactic centre as it logically can, namely at the cusped orbit, 
interior to which x\ orbits become self-intersecting. Consequently, 
a flow consistent with the BGSBU picture cannot be obtained with 
either a low sound speed or poor spatial resolution. It seems that 
previous simulations lacked the requisite resolution. We find that a 
consistent flow can be obtained for c s ~ lOkms -1 and grid spacing 
x < 5 pc. 

BGSBU did not provide a satisfactory explanation of the 
parallelogram-like structure of CO emission in the (/,v) plane. We 
find that the shocks form two sides of the CO parallelogram and 
conjecture that the prominence of the CO parallelogram is due 
to efficient conversion of atomic gas into molecular gas. Unfortu¬ 
nately, we do not follow the ISM’s chemistry. 

In our highest-resolution simulations the flow in the transi¬ 
tion region between the x\ and X 2 orbits is unsteady. We think this 
unsteadiness is probably a real physical phenomenon rather than 
a computational artifact, and is essentially turbulence generated in 
the region of high shear behind the shocks. We consider this un¬ 
steadiness, in conjunction with efficient conversion of atomic gas 
to molecular form in the shocks, provides a promising explanation 
of the observed asymmetry in the distribution of CO emission ei¬ 
ther side of the Galactic centre. 

While our simulations do provide strong support for the 
BGSBU picture, they do not explain all aspects of the observed HI 
and CO emission. There is, however, every prospect that further 
high-resolution simulations of flows in potentials similar to that 
used by BGSBU will explain all significant features. In this con¬ 
nection, two very worthwhile upgrades of our simulations would 
be an increase in the quadrupole moment of the bar, and inclusion 
of the conversion of gas between atomic and molecular forms. 
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